'''
Created on 08.06.2012

@author: Julia
'''

from __future__ import division
import numpy as np

c=[2]
print len(c)
class block():
    def __init__(self):
        self.H11=2.0
        self.L_inner=1.0

class system():
    def __init__(self):
        self.b1=block()
        self.b2=block()
        self.b3=block()
        self.b4=block()
        #self.psi=np.zeros((4), dtype=float)
    
    
    def get_gs(self):
        H=np.array([[self.b1.H11, -self.b1.L_inner, 0, 0], [-self.b1.L_inner, self.b2.H11, -1.0, 0], [0, -1.0, self.b3.H11, -self.b4.L_inner], [0, 0, -self.b4.L_inner, self.b4.H11]])
        #print H
        [u, v]=np.linalg.eig(H)
        idx=u.argsort()
        u=u[idx]
        v=v[:, idx]
        if v[:, 0].sum()<0:
            v[:, 0]=v[:, 0]*(-1)
        self.psi=v[:, 0]
        #print "psi:", self.psi
        self.energy=u[0]

"""
class wave_function():
    def __init__(self):
        self.v=np.zeros((4), dtype=float)
"""    
def new_basis(v, lr):
    if lr=="Left":
        a=v[0]
        b=v[1]
    else:
        a=v[2]
        b=v[3]
    norm=np.sqrt(a*a+b*b)
    a=a/norm
    b=b/norm
    return np.array([a, b])
        
    
        
def new_left(b1, b2, bas):
    res=block()
    res.H11=bas[0]*bas[0]*b1.H11+2*bas[1]*bas[1]-2*bas[0]*bas[1]*b1.L_inner
    res.L_inner=bas[1]
    return res

def new_right(b1, b2, bas):
    res=block()
    res.H11=bas[1]*bas[1]*b2.H11+2*bas[0]*bas[0]-2*bas[0]*bas[1]*b2.L_inner
    res.L_inner=bas[0]
    return res


def main():
    siteblock=block()
    length=raw_input("Length:")
    length=int(length)
    nsweeps=raw_input("Number of iterations:")
    nsweeps=int(nsweeps)
    allblocks=[]
    
    # Warmup sweep:
    allblocks.append(siteblock)
    for i in range(int(length/2)):
        S=system()
        S.b1=allblocks[i]
        S.b2=siteblock
        S.b3=siteblock
        S.b4=allblocks[i]
        S.get_gs()
        print i+1, S.psi[1], S.energy, 0
        basis=new_basis(S.psi, "Left")
        #bl=block()
        bl=new_left(allblocks[i], siteblock, basis)
        allblocks.append(bl)
    print "warmup sweep:", len(allblocks)
    #print int(length/2-1)
    bl=allblocks[int(length/2-2)]
    allblocks.append(siteblock)
    allblocks.append(bl)
    #allblocks.append(bl)
   
    #finite system sweep:
    for swp in range(nsweeps):
        # I discard the whole Reflect thingy as it is not important in this case and I dont understand it
        allblocks[int(length/2+1)]=allblocks[int(length/2-2)]
        for i in range(int(length/2+1), 2, -1):
            S=system()
            S.b1=allblocks[i-3]
            S.b2=siteblock
            S.b3=siteblock
            S.b4=allblocks[i]
            S.get_gs()
            print i-1, S.psi[2], S.energy, swp-0.5
            basis=new_basis(S.psi,"Right")
            bl=new_right(siteblock, allblocks[i], basis)
            allblocks[i-1]=bl
        #print i, "left sweep:", len(allblocks)
        print "blub", S.psi[0], S.energy, swp
        for i in range(0, int(length/2-1)):
            S=system()
            S.b1=allblocks[i]
            S.b2=siteblock
            S.b3=siteblock
            S.b4=allblocks[i+3]
            S.get_gs()
            print i+1, S.psi[1], S.energy, swp
            basis=new_basis(S.psi, "Left")
            bl=new_left(allblocks[i], siteblock, basis)
            allblocks[i+1]=bl
        print i, "right sweep:", len(allblocks)
    #energy=np.zeros((int(length)), dtype=float)
    #for i in range(50):
    #    energy[i]=allblocks[i].energy
    print "Length allblocks:", len(allblocks)
    return 0
#def get_groundstate(b1, b2, b3, b4, p):

main()
  
    